Simulating strongly-correlated multiparticle systems in a truncated Hilbert space 
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Representing a strongly interacting multiparticle wave function in a finite product basis leads to 
errors. Simple rescaling of the contact interaction can preserve the low-lying energy spectrum and 
long-wavelength structure of wave functions in one-dimensional systems and thus correct for the 
basis-set truncation error. The analytic form of the rescaling is found for a two-particle system 
where the rescaling is exact. A detailed comparison between finite Hilbert space calculations and 
exact results for up to five particles show that rescaling can significantly improve the accuracy of 
numerical calculations in various external potentials. In addition to ground-state energies, the low- 
lying excitation spectrum, density profile and correlation functions are studied. The results give a 
promising outlook for numerical simulations of trapped ultracold atoms. 
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I. INTRODUCTION 

Numerically solving the many-particle Schrodinger 
equation in one spatial dimension with short-range inter- 
action requires a truncated Hilbert space to make the cal- 
culations computationally tractable. This introduces in- 
accuracies especially in strongly correlated regimes, sys- 
tems with many degeneracies, and near avoided crossings. 
Here, we demonstrate a simple method that significantly 
improves results by rescaling the interaction strength. 

Understanding the complexities of many-body quan- 
tum physics remains a grand challenge, in particular with 
regard to excited states and dynamical problems. Ul- 
tracold atoms have recently become a testbed for both 
analytic and computational approaches to understanding 
quantum many-body physics. These systems can be pre- 
pared with few to millions of particles and the experimen- 
tal control over interactions, system size, and quantum 
states is rapidly improving [1-3]. 

However, exact theoretical results are rare, espe- 
cially for few-particle systems. Approximate approaches, 
like the mean-field Gross-Pitaevskii equation, accurately 
model a large number of weakly interacting atoms. Now, 
there is a considerable amount of research that focuses on 
the creation of entanglement in ultracold atomic systems. 
These require calculations at a single-particle level. Much 
success with fixed basis-set expansions has been achieved 
using Wannier functions of an optical lattice leading to 
the Bose-Hubbard model [4] and with harmonic oscillator 
eigenfunctions [5] or a plane-wave basis [6, 7] for trapped 
systems. These approaches are limited by accessing only 
a finite Hilbert space. Methods like time-evolving block 
decimation [8] and time-dependent density matrix renor- 
malization [9, 10] exploit spatial entanglement properties 
to allow for time-dependent simulations of larger multi- 
particle problems. 

A different approach is taken by the multiconfigura- 
tional time-dependent Hartree (MCTDH) method [11- 



13]. It variationally optimizes a set of single-particle 
functions, that define the accessible Hilbert space at any 
given moment in time. It has been successfully applied 
to ultracold gases [14-16]. Even so, the truncated ba- 
sis leads to errors in strongly interacting systems. Here, 
we analyze the origin of errors that arise from basis-set 
truncation, suggest a scheme to compensate for them, 
and benchmark the performance of the scheme. 

Trapped bosonic atoms confined to a single spatial di- 
mension undergo a crossover from a Bose-Einstein con- 
densate to the fermionized and strongly correlated Tonks- 
Girardeau (TG) gas [17] as repulsive short-range interac- 
tions are tuned from weak to very strong. This system 
is well suited to benchmark computational many-body 
methods as exact results are available from integrable 
models in the limits of infinite interactions (TG [17]) 
or vanishing external potential (Lieb-Liniger model [18]). 
Previous calculations of the ground state with exact di- 
agonalization in a fixed product basis [5, 19, 20] and with 
the MCTDH method [14, 21] have shown a peculiar fea- 
ture: Specific finite values of the interaction strength in 
these calculations provide a reasonable approximation of 
the completely fermionized state, which in exact theory 
is only reached for infinite interaction. Our calculations 
show that larger values of the interaction produce worse 
results in the approximate calculation for the energy and 
lead to spuriously enhanced density oscillations. In this 
paper we provide a clear explanation of this effect based 
on exact results for two particles. 

We find that the origin of the problem lies in the 
truncation of the Hilbert space. Previously, elaborate 
schemes have been developed to generate effective Hamil- 
tonians that remove this problem, however these tech- 
niques are cumbersome and require careful application 
[22-25]. A simpler, yet powerful method is desirable to 
accurately describe different interacting regimes for a va- 
riety of applications [26-29] . Here we present a rescaling 
method that increases accuracy or decreases the required 



2 



Hilbcrt space by mapping the interaction in the simula- 
tions to the physical interaction. Therefore numerical 
simulations can be sped up and used to approach previ- 
ously unaccessible regimes and systems. 

In this paper we initially present an exact rescaling 
procedure for two particles in a one-dimensional peri- 
odic box as described previously in Ref. [7]. A more 
general scheme is proposed that only requires the knowl- 
edge of the numerical value of the interaction strength at 
the TG energy. This approach is closely related to the 
known renormalization of the two-body T matrix but ex- 
tends this method into the strongly interacting regime. 
The approach is tested for different system geometries, 
particle numbers and external potential shapes. The re- 
sults are compared to analytical results for two atoms in 
a harmonic potential [30] , the Lieb-Liniger description of 
atoms confined to a one-dimc!nsional periodic system [18] 
and infinitely strong interacting atoms in one dimension 
described by a TG gas [17]. 

II. THEORETICAL BACKGROUND 

A. Finite basis-set expansion 

Consider N interacting quantum particles in one di- 
mension subject to an external potential \4xt(a^) t), which 
may depend on time: 

N 

Hmb ^^hi{xi) + ^gd{xi- Xj), (1) 

i—l i<j 

where hi = —■^-^+Vcxt{xi,t) is the one-body Hamilto- 

nian and for simplicity we assume that all particles have 
the same mass m and interact by contact interactions of 
strength g. This model describes, e.g., quantum gases 
of ultracold atoms in an elongated and tightly confin- 
ing trapping potential [31]. Although we specifically deal 
with identical bosons later on, at this point considera- 
tions are not restricted to a specific quantum statistics 
and the Hamiltonian (1) could describe single- or multi- 
component Bose or Fermi gases or mixtures thereof. All 
these possibilities are of interest and related to actual or 
possible experimental scenarios. 

The true A^-particle wave function can be expanded in 
the form 

'^{xi,...,XN,t) = ^Aj(i)$j(xi,...,Xjv), (2) 

J 

in a basis consisting of products of single-particle wave 
functions 

N 

^j{xi,...,Xn) = Yl(l)jk{Xk), (3) 

fe=l 

where we have introduced the multi-index J = 
(ji, . . . , jat) for the set of N single-particle indices. The 



single-particle wave functions (or mode functions) are 
mutually orthonormal by / (j)*{x)(pk{x)dx = 5jk- 

In practice we chose a finite set of M single-particle 
functions to define the finite basis-set expansion. The size 
of the finite basis of multiparticle Hilbcrt space is nom- 
inally , although for identical particles this number 
can be significantly reduced by accounting for bosonic or 
fermionic exchange symmetry in a well-known manner 
by symmetrizing the products to permanents or Slater 
determinants, respectively. 

The standard procedure for approximating the full 
problem (1) is to simply truncate the Hilbcrt space to 
span a finite basis and represent the Hamiltonian as a 
finite matrix with elements 

Hjj = ($/|i/,„fa|$j). (4) 

We denote the solution of the corresponding matrix 
Schrodingcr equation for the coefficient vector with Aj 

and the eigenvalues of the truncated matrix E^,. These 
approximate the exact coefficient vector Aj and eigenval- 
ues E^, which form part of the spectrum of Hmb- Since 
the basis-set truncation can be understood as a varia- 
tional procedure with a restricted variational space, the 
approximate ground-state energy Eq > Eq is an upper 
bound for the true one. By increasing the number of 
modes M and thereby computational space, the approx- 
imation improves and the approximate energy converges 
to the correct one. However, for an interacting system 
this convergence is painfully slow, as seen in the inset in 
Fig. 1. 

The purpose of this paper is to improve the approx- 
imations for a given, finite, basis-set expansion. While 
previous works have explored elaborate procedures to ap- 
proximate an effective finite Hamiltonian with the exact 
eigenvalues Ei, [22-24], we show that it is possible to 
obtain significant improvements by simply rescaling the 
interaction parameter g. 



B. Two particles on a ring 

For two particles without external potential an analytic 
expression for the rescaled interaction constant g can be 
found that reproduces the correct ground-state energy 
and the Fourier components of the wave function up to 
a cutoff. The rescaling procedure is thus exact for this 
case. 

Consider N = 2 bosons in a one-dimensional box with 
length L and periodic boundary conditions. This model 

might be realized with ultracold atoms in a tightly fo- 
cused ring trap. The Hamiltonian (1) simplifies to 

^^- = -2;;ite + a^j+^^("^-"^) 

where Xi is the spatial coordinate of atom i and g the 
physical interaction strength. A finite basis-set expansion 
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FIG. 1. Relative deviation of the approximate ground-state 
energy from the exact Lieb-Liniger (LL) result [18] before 
rescaling (white area of panel a) and after rescaling (gray 
area of panels a and b) for five particles in a box of length 
L with periodic boundary conditions. Rescaling reduces the 
maximum errors from 23% to 0.6% (M — 13). The inset 
shows the relative energy deviation vs the number of modes 
M for three particles ai g = 1.9099. Results for the trun- 
cated Hamiltonian without rescaling (squares) converge very 
slowly while rescaling (circles) significantly reduces the errors 
already for small M. 



for this system is 



(M-l)/2 

E 

fci,fc2 = -(A/-l)/2 



-kixi 



-k2X2 



(6) 



where we assume the number of momentum modes M in 
the expansion to be odd, for simplicity. First consider 
the case for the full Hilbert space with M = oo. To 
find the ground state we substitute Eq. (6) and Eq. (5) 
into the Schrodinger equation {E — H2p)ip — 0. Pro- 
jecting to zero-momentum solutions by multiplying with 
exp[— i27rfc(a;i — X2)/L] and integrating over the particle 
coordinates leads to 



{E - 2Eok^)Ck,-k 
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L 



E 



Cq-q. 



(7) 



Here, £^o = 27r^?i^/(mL^) is the smallest nonzero single- 
particle energy. The right-hand side of Eq. (7) is inde- 
pendent of k and therefore C-q^q oc {E — 2Eoq^)~^ and 
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(8) 



which relates the exact energy E to the physical interac- 
tion strength g. Employing the finite basis-set expansion 



limits the sum over q to only M terms and thus yields a 
different, approximate, value E for the energy. Here, we 
choose a different path and introduce the rescaled inter- 
action strength g by 
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2q2 



(9) 



(10) 
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where the last equation defines the constant go- While 
Eq. (9) guarantees that g yields the exact energy E it is 
Eq. (10) that provides the rescaled interaction constant 
5 as a function of g: 

9 
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1 + g/90 ' 



(11) 



Therefore, solving the Schrodinger equation with this 
rescaled value g instead of g in the Hamiltonian (5) in 
the finite basis-set expansion gives the exact energy E 
for the physical interaction strength g. Also the finite 
number of expansion coefficients Ck.~k is identical to the 
full Hilbert space expansion and, according to Eq. (7), is 
given by 

Ck.-k = ^ f ^ , , , (12) 



E - 2Eok^ ' 



where 



A = i 



q— — oo 



^1,-q ^ ^ E ^1-1- 



(M-l)/2 



(13) 



+ 0(Af-5). 



(14) 



The value of g^ is found by expanding the sums in powers 
of 1/M from 

EqL 2 2E-Eq 
~ M + 2EoAP 

For a large number of modes and sufficiently low energy 
E [with the approximate condition E/{2Eo)] we 

thus find to leading order that go is independent of the 
energy and given by 
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-MEqL. 
2 ' 



(15) 



Figure 2 shows how the rescaling in approximation (15) 
significantly improves the ground-state energy. For AI — 
6 the energy after rescaling is only 3% off the exact energy 
whereas before rescaling the difference was ^ 20%. 

The constant go is readily interpreted from Eq. (11) as 
the value of the rescaled interaction g where the physi- 
cal interaction g reaches infinity. Physically, this is the 
TG limit where the contact interactions are arbitrarily 
strong and the system is strongly correlated. Our results 
thus indicate that a finite basis-set expansion with a fi- 
nite value of the rescaled interaction constant g = go can 
produce exact results for infinite values of the physical in- 
teraction strength. Since the energies and wave functions 
in the TG limit can also be calculated by the Bose-Fermi 
mapping [17], this limit is an important reference point. 
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FIG. 2. Ground-state energies for two interacting bosons in 
a periodic box before (lines) and after (symbols) rescaling ac- 
cording to Eqs. (11) and (15) compared with the exact (Lieb- 
Liniger) energy (solid line) . The calculated energies have been 
obtained via exact diagonalization of the truncated Hamilto- 
nian matrix with different numbers of modes M. Clearly, 
rescaling of the interaction reduces the errors in the energies 
significantly. 



C. Empirical rescaling 

The analytical results of the preceding section moti- 
vate us to extend the idea of rescaling to more than two 
particles and to problems with arbitrary external poten- 
tial, which, generally, cannot be solved analytically. We 
propose to replace the Hamiltonian (1) by the rescaled 
version 



N 



Hmb = ^hi{xi) + ^g5{xi - Xj), 



(16) 



where we have only changed the value of the interac- 
tion strength from the physical value g. Within a fi- 
nite basis-set expansion we obtain a matrix with elements 
Hjj, which is to be used for numerical simulation. The 
smallest eigenvalue is the approximate ground-state en- 
ergy Ecig)- For the rescaled interaction g we continue to 
use Eq. (11) and determine the value of go by requiring 
that the approximate ground-state energy at go equals 
the exact energy of the TG limit: 



Eg {go) = Etg- 



(17) 



A finite solution for go can always be found since the 
approximate value of Eq {g) overestimates the real value 
and Eg {g) is a monotonously growing function of g with 
Eg{0) < -Etg- The value of Etg is equal to the ground- 
state energy of a system of noninteracting fermions due to 
the Bose- Fermi mapping theorem [17]. It is given by the 
sum of single-particle energies, which can be determined 
analytically or numerically to very high precision. 



This rescaling procedure is very simple to implement, 
since only knowledge of the value of the single constant 
go is needed. We assume here for simplicity that go is 
independent of the energy of the wave function and can 
be used for any set of (low-lying) excited states as well 
as time-dependent processes involving such states. As 
an approximation, this is consistent with the analytic 
findings for two particles from the previous section. 



D. Connection to T matrix renormalization 

The rescaling approach introduced in Sec. II C is 
closely connected to the well-known renormalization of 
the scattering T matrix (for a readable account see 
Ref. [32]). In this approach, the coupling constant of an 
effective many-body theory is renormalized to produce 
the correct T matrix for the two-body scattering prob- 
lem, which is known exactly. Indeed, Eq. (11) for the 
rescaled interaction strength together with Eq. (15) for 
the value of go are identical to the result of Ref. [32] for 
the renormalized T matrix of a high-momentum cutoff in- 
troduced there due to discretization of space (identifying 
the computational volume of momentum space 2'jtM/L 
with the Brillouin zone of Ref. [32]). This result was 
rigorously derived for weak interactions and high cut- 
offs. In Sec. II B we have extended the approach to the 
strongly interacting regime for two particles where now 
the rescaled interaction constant in principle becomes 
energy dependent. The extension to nonperturbative, 
strongly interacting multiparticle systems in Sec. II C 
with Eq. (8) for the value of go is heuristic in nature. The 
purpose of the numerical studies reported in the follow- 
ing section is to establish the usefulness and to quantify 
the accuracy of this approach. 

Rescaling improves numerical calculations by either in- 
creased acciiracy or decreased computational effort for a 
given desired accuracy (CPU time ~ for our five par- 
ticle simulations). For g > go we find unphysical results. 
Therefore, CPU times are further reduced as only results 
for g < go arc used in the rescaling. 

Before detailing our numerical results and benchmark- 
ing, we introduce the MCTDH method suitable for gen- 
eral trapping potentials. 



E. Variationally optimized basis 

The finite basis-set expansion (2) and (3) expands the 
multiparticle wave function in products of a given or- 
thonormal set of single-particle functions. Optimizing 
this set variationally is the aim of the MCTDH method. 
This is achieved with a time-dependent variational ansatz 
with a fixed number of basis functions without further 
approximation [11], in contrast to the density matrix 
renormalization group approach, which further decom- 
poses the problem into bipartite systems [9, 10] before 
finding an optimized basis. 
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The single-particle wave functions (f>j{x,t) arc now al- 
lowed to be time dependent in addition to the time de- 
pendence of the coefficients Aj{t). The variational equa- 
tions of motion for the coefficients and the single-particle 
wave functions are discretized in time (and space) and 
solved numerically with the MCTDH program package 
[12, 13]. 

Eigenenergies and eigenstates can be obtained by the 
relaxation method [33], i. e., by propagation in negative 
imaginary time. Here we use a modification, improved re- 
laxation [34] , which can also be used in block form [35] to 
simultaneously compute a set of eigenvectors. For a fixed 
basis of single-particle functions, the MCTDH method 
reduces to the exact diagonalization method. There- 
fore, the key advantage is that the MCTDH method 
always stays in a variationally optimized basis. This 
generally ensures higher accuracy or keeps the Hilbert 
space as small as possible for a given accuracy. The 
MCTDH method is therefore effective at describing dif- 
ferent trapping potentials. However, it requires a numer- 
ically unattainable number of single-particle functions to 
describe strongly interacting systems. 

In practice, the single-particle wave functions are rep- 
resented on a large number of primitive basis functions, 
the choice of which depends on the spatial geometry (typ- 
ically 128 and 162 grid points in the harmonic oscilla- 
tor discrete variable representation for the harmonic and 
the double-well potentials respectively and 155 points in 
the exponential discrete variable representation as per 
Ref. [12]). The setup of the MCTDH method for ultra- 
cold atom calculations follows Ref. [14]. In particular, 
the contact potential is replaced by a narrow Gaussian 
for practical purposes. We have carefully studied the 
dependence of our calculations on the parameters of the 
primitive basis and the Gaussian interaction to make sure 
that the results presented here are converged and the in- 
fluence of these further approximations is negligible. 



III. BOSONS IN A RING 

We first benchmark the rescaling method of Sec. II. C. 

for the Lieb-Liniger model [18] of bosons in a box with 
periodic boundary conditions, where exact eigenstate en- 
ergies are easily available for any values of the interaction 
strength and particle number. Figure 1 shows the relative 
deviation from the exact ground-state energy obtained 
from MCTDH calculations before and after rescaling. 
Rescaling is seen to significantly reduce the error. The 
rescaled results reproduce the exact energy in the TG 
limit by construction but the improved accuracy for fi- 
nite interaction strength is nontrivial. As expected, both 
the raw and rescaled data improve with increasing M. 
The maximal errors for, e.g., M = 13 can be reduced 
from ~ 23% down to ~ 0.6%, i. e. a factor of 40. This is 
a significant improvement. 

In order to test how well the wave function is approx- 
imated in the TG limit we investigate the reduced pair 
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FIG. 3. Density-density correlation function for five particles 
and M = 13 single-particle functions in a periodic box. This 

plot compares the reduced two-particle density pif^ of Eq. (18) 
for different interaction strengths from MCTDH calculations 
without rescaling with the exact TG expression (the area be- 
low this curve is shaded). We find a good agreement with the 
exact solution close to g'o = 5.3958EoL. This is shown in the 
inset where the overlap F^<^\p^^\ pl'^) of Eq. (19) is plotted 
as a function of the interaction g. There, the maximum occurs 
close to go. 

density 

p^i\Q,x)= j-- J dX3...dXN\'^<^a\0,x,X3,...,XN,t)\^■ 
{^S) 

The results of MCTDH calculations are shown in Fig. 3 
for different values of the interaction g along with the 
exact solution in the TG limit /3|^'^(0,:r) = [(2 -|- 
cosa;) sin^ |]/37r^ [17]. Agreement between the finite 
basis-set expansion at the value of go = 5.3958£'o^ and 
the exact solution is reasonable. The discrepancies can 
well be explained by the absence of short-wavelength 
modes in the finite expansion. Most importantly, the 
agreement between the approximate and exact functions 
clearly worsens for interaction strengths larger than go, 
where the structure in the finite basis-set calculation 
overshoots. 

In order to quantify the quality of the approximation, 
we define the overlap with the exact solution at the TG 
regime by 

Jp^^i0,x)pi'\0,x)dx 

^J J pf (0, x)pf (0, x)dx^\ p^^'iO, x)p™(o, ^)dx ' 

(19) 
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FIG. 4. Maximal relative deviation of the energy after rescal- 
ing from the exact result Ell over the range of interaction 
strengths < g < oo. The particle number and box length L 
are varied while the density is kept fixed. 
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FIG. 5. Density of five particles in a harmonic potential for 
different interaction strengths without rescaling and AI — 13. 
Analogous to Fig. 3 the best agreement with the exact TG 
result, which is given by the shaded areas, can be found at 
g = go, where go = 15.516EoL. 



which is shown in the inset of Fig. 3. We find that the 
maximum of F^^^ is reached near go, which justifies the 
choice of go through Eq. (17). 

In order to study the scaling with particle number, 
we have performed exact diagonalization calculations for 

= 2 . . . 6 particles in a ring of length L = N/h, where 
the particle density h was fixed to the same value as in 
Fig. 3. Figure 4 shows the maximal relative deviation 
of the rescaled energy as a function of particle number. 
The results show that retaining M = 27V — 1 single- 
particle functions gives a maximal relative error bounded 
by \E/Ell - 1| < 1-5% for any value of g > 0. Ahhough 
for clarity only odd numbers of single-particle modes M 
are shown in Fig 4, we have verified that these findings 
also hold for even M. 

If we assume that the scaling M = 2N — 1 is sufficient 
to maintain accuracy uniformly over the whole range of 
interaction strengths also for larger particle numbers, we 
still find that the numerical effort of MCTDH or exact 
diagonalization studies grows exponentially with N. This 
limits studies of this sort to fewer than 10 particles with 
current-day computers. 



IV. BOSONS IN A HARMONIC POTENTIAL 

In order to test the rescaling method in the presence 
of an external potential we consider a harmonic trapping 
potential 



(20) 



ing Eo,H = fiw and Lh — y respectively, where w is 
the harmonic oscillator frequency. 

Figure 5 compares the single-particle density of the TG 
wave function with MCTDH simulations for five parti- 
cles at different values of the interaction strength with- 
out rescaling. The agreement between the exact and 
the finite basis-set results at the value of the interaction 
strength go = 15.516 . . .Eq^hLh is much better than for 
larger values of g. This result is analogous to the ring 
case and supports our proposition that go from Eq. (17) 
provides a good approximation for the wave function. 

Full analytical solutions are available for two parti- 
cles in a harmonic potential. The exact energies for the 
ground and excited states are found from the implicit 
equation [30] 
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Eo.hLh 



= 2V2' 



2Eo 



+ 1) 



2Eo. 



(21) 



where F is the gamma function. Figure 6 compares 
the ground-state and a few excited-state energies with 
rescaled MCTDH results. The highest few states and the 
ground state have been chosen as they are more likely to 
show the largest errors within the block-diagonalization 
scheme (see below). We emphasize that a single value of 
go = 29.75 obtained from Eq. (17) was used for rescal- 
ing. The procedure not only maintains the correct or- 
der of excited states but also quantitatively describes the 
excited-state energies remarkably well. 



BOSONS IN A DOUBLE WELL 



For this and the following section we use as natural en- 
ergy and length scales the harmonic oscillator level spac- 



The rescaling method is further tested for the descrip- 
tion of excited states of the TG gas in a double-well po- 
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FIG. 6. Energies of the ground state and the sixth to eighth 
excited states with M = 15 of two bosons in a harmonic well 
after rescaling compared to the exact results from Eq. (21). 
The results clearly show the high accuracy of the rescaled 
results. The largest error occurs for the seventh excited state 
where it is decreased from ~ 2% before rescaling (not shown 
on this graph) to ~ 0.3% after rescaling. 

tential. Exact solutions are not available for interacting 
multiparticle systems in the double well except for the 
TG gas, where their exact energies are simply found from 
sums of single-particle energies (which can be determined 
with great accuracy) due to the Bose-Fermi mapping of 
Rcf. [17]. 

We consider the double- well potential 

V..t{x)/Eo,H = I + ^^"^^^/^"^^ (22) 

with varying barrier height h at the center. Figure 7 
shows the energies for the ground and excited states from 
rescaled MCTDH calculations and analytical results. All 
MCTDH results were obtained in the TG regime at g = 
go- The value of go has been determined for each value 
of h from a dedicated ground-state calculation. 

The results reported in Fig. 7 were calculated using 
the block-improved-relaxation method [35] to efficiently 
determine several states simultaneously. Since an iden- 
tical set of single-particle functions is used for the finite 
basis-set expansion for all states, the result is not as well 
variationally optimized for each individual state as the 
single-state calculation that was used to determine go- 
This explains the discrepancy of the numerical and exact 
ground-state energies seen in Fig. 7, which also provides 
a convenient order-of-magnitude estimate for the error 
of the block-improvcd-rclaxation scheme. A shift to high 
energies for the whole set of numerical values compared 
to the exact ones of the same approximate magnitude is 
clearly observed in the data. 

In addition, the second and third excited states in 
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FIG. 7. Excited states of the Tonks- Girardeau (TG) gas in 
a double-well trap. Shown are the energies of the ground 
state and the first four excited states minus the exact TG 
ground-state energy for five bosons in a double-well poten- 
tial vs barrier height h/Eo^H- The graph compares rescaled 
MCTDH calculations with M = 13 single-particle functions 
(solid line with bullets) and exact TG energies (dashed line). 

Fig. 7 should be degenerate at ft- = 0, which is not re- 
produced by the rescaled MCTDH results. By the Bose- 
Fermi mapping, these two states are identified as configu- 
rations (1, 1, 1, 1, 0, 0, 1,0,.. .) and (1, 1, 1, 0, 1, 1, 0,0,.. .) 
in the occupation number basis of fermions in the har- 
monic oscillator {h = 0) potential. The two states in- 
volve different modes in the fcrmionic description and 
thus have different sensitivity to the finite set of available 
single-particle functions in the block-improved-rclaxation 
method. This is a problem that the proposed rescaling 
method cannot fix. 

The problem vanishes when more single-particle func- 
tions are used, but this is costly. Single-state simula- 
tions indeed show that the energy difference approaches 
zero for an increasing number of single-particle functions. 
These problems do not affect the validity of rescaling 
which is shown to work adequately and to be independent 
of the shape of the external potential. The energy dif- 
ferences compared to the analytical results range around 
about 1%. 

For h = lOEo^H a near degeneracy between the 
third and fourth state emerges. In contrast to the 
previously discussed degeneracy at vanishing barrier, it 
is very well described by the simulation. For large 
enough barrier height the double well resembles two 
separated anharmonic wells. According to the Bose- 
Fermi mapping, we identify the degenerate states by 
the superpositions of (1, 1, 1,0, . . .); (1,0, 1,0, . . .) and 
(1, 0, 1,0,.. .); (1, 1, 1,0,.. .), showing the fermionic con- 
figurations in the left- and right-hand wells, respectively. 
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Since these two states involve similar fermionic modes, 
their degeneracy is well described by the truncated basis- 
set expansion with a limited number of single-particle 
functions. 

Overall, we find good description of the excited states 
by the roscaled simulation in the double well for a large 
variation of barrier height h from a harmonic trap {h = 0) 
to almost completely separated wells. 

VI. CONCLUSION 

We have proposed a simple rescaling procedure for 
many-body calculations with a contact interaction in 
one dimension in a truncated Hilbert space based on 
Eqs. (11) and (17). The method requires knowledge 
of only a single parameter. This is found from the inter- 
action strength where the numerical ground-state energy 
matches the TG energy, which is easily found from single- 
particle calculations. For two particles in a ring geometry 
the rescaling exactly reproduces the correct energies and 
expansion coefficients of the wave functions. Empirical 



evidence suggests that the scheme can be extended to 
multiple particles, external trapping potentials, and low- 
lying excited states. Testing the method against exact 
results for ring geometries and harmonic and double-well 
potentials shows significant improvements in accuracy. 

Because rescaling works for excited-state calculations 
and for arbitrary potentials it is expected to significantly 
improve time-dependent simulations as well, as long as 
the state of the system is well described by the low-lying 
part of the energy spectrum. 
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